#### "Economic Inequality, Income and political participation in authoritarian regimes" ####
# authors: "Pelke, Lars"
# date: 2020-06-08
# written under "R version 3.6.0 (2019-07-05)"

#### Preliminaries ####

R.version$version.string

# clear workspace
rm(list=ls())

# set working directory

# loading packages

library(countrycode)
library(tidyverse)
library(viridis)
library(haven)
library(lme4)
library(sjstats)
library(texreg)
library(interplot)
library(ggrepel)
library(stargazer)
library(ggpubr)
library(dotwhisker)
library(broom.mixed)


#### Import Data ####
wvsdata <- readRDS("data/wvsdata.rds") 

#### Reducing to authoritarian regimes ####

wvsdata <- wvsdata %>%
  filter(v2x_regime_amb<=4)

#### Analyze Regime with elections multiparty:
# at least "constrained. at least one real opposition party is 
# allowed to context but competition is highly constrained -legally or informally. 

wvsdata <- wvsdata %>%
  filter(v2elmulpar_ord >=2)

wvsdata <- wvsdata %>%
  mutate(country_year = stringr::str_c(country_name, year, sep=" "))

table(wvsdata$country_year, is.na(wvsdata$gini_mkt))
table(wvsdata$country_year, is.na(wvsdata$v2elvotbuy))

number_observations <- wvsdata %>%
  count(country_year)

#### Building variables ####

#Replacing missings in Dependent Variables

wvsdata <- wvsdata %>%
  filter(E264>=0 | E257>=0) %>%
  mutate(vote = ifelse(E264==1, 1,
                       ifelse(E264==2, 1,
                              ifelse(E257==1, 1, 0))))
table(wvsdata$vote)
summary(wvsdata$vote)

number_observations <- wvsdata %>%
  count(country_year)

#Replacing missings in independent variables

wvsdata <- wvsdata %>%
  filter(X001>=0)%>%
  mutate(sex= X001)

wvsdata <- wvsdata %>%
  filter(X003>=0)%>%
  mutate(age= X003)

wvsdata <- wvsdata %>%
  filter(X025>=0)%>%
  mutate(education= X025)

wvsdata <- wvsdata %>%
  filter(X028>=0)%>%
  mutate(employement_status= X028)

wvsdata <- wvsdata %>%
  filter(X047>=0)%>%
  mutate(income_deciles= X047)

wvsdata <- wvsdata %>%
  mutate(income_quintiles = ifelse(income_deciles==1, 1,
                                    ifelse(income_deciles==2, 1,
                                           ifelse(income_deciles==3, 2,
                                                  ifelse(income_deciles==4, 2,
                                                         ifelse(income_deciles==5, 3,
                                                                ifelse(income_deciles==6, 3,
                                                                       ifelse(income_deciles==7 , 4,
                                                                              ifelse(income_deciles==8, 4, 5)))))))))


wvsdata <- wvsdata %>%
  filter(X011>=0)%>%
  mutate(children= ifelse(X011>0, 1, 0))

wvsdata <- wvsdata %>%
  filter(X028>=0)%>%
  mutate(unemployed= ifelse(X028==7, 1, 0))

table(wvsdata$vote)
# 9691 non vote, 29915 vote

number_observations2 <- wvsdata %>%
  count(country_year)

## Inspect which countries are delete because missing obserations at macro-level variables ##

setdiff(number_observations$country_year, number_observations2$country_year)

###Data Manipulation of GINI measures, fill some missing value with last national observations

missing <- wvsdata %>% 
  group_by(country_year) %>% 
  summarise(sumNA = sum(is.na(gini_mkt)))
missing

wvsdata <- wvsdata %>%
  mutate(gini_mkt=ifelse(country_year=="Algeria 2014", 36.7, gini_mkt)) %>% # Algeria 2011
  mutate(gini_disp=ifelse(country_year=="Algeria 2014", 35, gini_disp)) %>%
  mutate(gini_mkt=ifelse(country_year=="Azerbaijan 2011", 38.3, gini_mkt)) %>% # Azerbaijan 2011
  mutate(gini_disp=ifelse(country_year=="Azerbaijan 2011", 30.3, gini_disp)) %>%
  mutate(gini_disp=ifelse(country_year=="Lebanon 2013", 38.5, gini_mkt)) %>% # Lebanon 2012
  mutate(gini_disp=ifelse(country_year=="Lebanon 2013", 36.5, gini_disp))

wvsdata <- wvsdata %>%
  mutate(gini_mkt=ifelse(country_year=="Haiti 2016", 57.6, gini_mkt)) %>%
  mutate(gini_disp=ifelse(country_year=="Haiti 2016", 53.7, gini_disp)) %>%
  mutate(gini_mkt=ifelse(country_year=="Nigeria 2012", 46.8, gini_mkt)) %>%
  mutate(gini_disp=ifelse(country_year=="Nigeria 2012", 44.1, gini_disp)) 

wvsdata <- wvsdata %>%
  drop_na(gini_mkt)
rm(missing)

#### Previous democratic experience ####

table(wvsdata$democratic_expericene) 
table(wvsdata$democratic_expericene_all)

democratic_exp <- wvsdata %>%
  filter(democratic_expericene_bi ==1) %>%
  group_by(country_name) %>%
  count()
# authoritarian regime that has previous democratic experience: Armenia (1990-1994), Belarus (1991-1996), Palestine (2004-2006), Russia (1992, 1993), 
# Thailand (1998-2005, 2012)

rm(democratic_exp)

#### Social Group Base Regime ####

table(wvsdata$regime_support_group_max) 

#### Country and Country-Year overview ####

country_year_observations <- wvsdata %>% 
  group_by(country_year) %>% 
  summarise(sumNA = sum(is.na(gini_mkt)))
country_year_observations
rm(country_year_observations)

country_observations <- wvsdata %>% 
  group_by(country_name) %>% 
  summarise(sumNA = sum(is.na(gini_mkt)))
country_observations
rm(country_observations)

count <- wvsdata %>%
  group_by(S025, income_quintiles) %>%
  tally() 
rm(count)


number_observations3 <- wvsdata %>%
  count(country_year)
setdiff(number_observations$country_year, number_observations3$country_year)

##################################
#Group and grand mean centering
##################################

# Author: Zoltan Fazekas
centFUN <- function(x) {
  x - mean(x, na.rm = TRUE)
}

# Median Cenetering for factor variables
centFUN_median <- function(x) {
  x - median(x, na.rm = TRUE)
}

#### Data manipulation ####

wvsdata <- wvsdata %>%
  mutate(democratic_expericene_all = ifelse(is.na(democratic_expericene_all), 0, democratic_expericene_all), 
         democratic_expericene = ifelse(is.na(democratic_expericene), 0, democratic_expericene))

# Use the function to grand-mean center the macro-level predictor
wvsdata$gini_mktCGM <- centFUN(wvsdata$gini_mkt)
wvsdata$gini_dispCGM <- centFUN(wvsdata$gini_disp)
wvsdata$v2xel_frefairCGM <- centFUN(wvsdata$v2xel_frefair)
wvsdata$e_migdppclnCGM <- centFUN(wvsdata$e_migdppcln)
wvsdata$v2elvotbuyCGM <- centFUN(wvsdata$v2elvotbuy)
wvsdata$regime_support_group_maxCGM <- centFUN(wvsdata$regime_support_group_max)
wvsdata$democratic_expericene_allCGM <- centFUN(wvsdata$democratic_expericene_all)
wvsdata$democratic_expericeneCGM <- centFUN(wvsdata$democratic_expericene)


summary(wvsdata$e_migdppclnCGM)
summary(wvsdata$v2xel_frefairCGM)
summary(wvsdata$gini_mktCGM)
summary(wvsdata$gini_dispCGM)
summary(wvsdata$v2elvotbuyCGM)
summary(wvsdata$regime_support_group_maxCGM)
summary(wvsdata$democratic_expericeneCGM)
summary(wvsdata$democratic_expericene_allCGM)



#Use the fuction to group-mean center the individual predictors
wvsdata <- wvsdata %>%
  group_by(S025) %>%
  mutate(ageCWC = centFUN(age))
summary(wvsdata$ageCWC)

wvsdata <- wvsdata %>%
  group_by(S025) %>%
  mutate(educationCWC = centFUN(as.numeric(education)),
         income_quintilesCWC = centFUN(as.numeric(income_quintiles)),
         income_decilesCWC = centFUN(as.numeric(income_deciles)))

summary(wvsdata$educationCWC)
summary(wvsdata$income_quintilesCWC)
summary(wvsdata$income_decilesCWC)

#Rescaling predicator variables
summary(wvsdata$educationCWC)
summary(wvsdata$sex)
summary(wvsdata$ageCWC)
summary(wvsdata$income_quintilesCWC)
summary(wvsdata$S003)
wvsdata$ageCWC <- wvsdata$ageCWC/10


##################################
#Descriptive Statistics
##################################

count <- wvsdata %>%
  group_by(S025, income_quintiles) %>%
  tally() 

wvsdata <- wvsdata %>%
  group_by(S025, income_quintiles) %>%
  mutate(vote_byincome = mean(vote))

wvsdata <- merge.data.frame(wvsdata, count)

rm(count)

## Figure B1 ##

ggplot(data=wvsdata, aes(x=income_quintiles, y = vote_byincome/n))+
  geom_col() +
  facet_wrap(~country_year, ncol = 5) +
  theme_bw() +
  xlab("Income Quintiles") +
  ylab("Voting (Yes)") +
  scale_y_continuous(labels = scales::percent) +
  labs(y = "Average Participation in authoritarian elections by income groups",
                                                     x = "Income Quintiles",
                                                     title = "",
                                                     caption = "Data from World Value Survey and V-Dem dataset, version 10")

ggsave("output/B1.pdf", width = 22 , height = 25 , units = c("cm"))

#### Min and Max Voting ####

count_votes <- wvsdata %>%
  group_by(country_name) %>%
  summarize(vote = mean (vote))

#### Scatterplot Inequality and Voting ####

wvsdata_country <- wvsdata %>%
  group_by(country_year) %>%
  summarize(voting_mean = mean(vote),
            inequality = mean(gini_mkt))

ggplot(data= wvsdata_country, aes(x=inequality, y = voting_mean, label = country_year)) +
  geom_point() +
  geom_smooth(method = lm, colour = "black") + 
  geom_text_repel(size = 3) +
  theme_classic() +
  labs(y = "Average Vote Share in Country Year",
       x = "Market Inequality",
       title = "",
       caption = "Data from World Value Survey and V-Dem dataset, version 10")

ggsave("output/A1_scatter_vote_inequality.pdf")

rm(wvsdata_country)

#### Histogramm Voting ####

ggplot(data=wvsdata, aes(x = vote)) +
  geom_histogram(binwidth = 1, fill = "white", colour = "black") +
  theme_classic() +
  scale_x_continuous(
    breaks = c(0, 1),
    labels = c("No", "Yes")) +
  labs(y = "Number of individuals",
       x = "Vote (yes/no)",
       title = "",
       caption = "")
ggsave("output/A1_histogramm_vote.pdf")

##Summary statistics ##

cols <- c("sex", "ageCWC", "educationCWC", "children", "unemployed", "income_quintilesCWC",
          "gini_mktCGM", "gini_dispCGM", "v2xel_frefairCGM", "v2elvotbuyCGM", "e_migdppclnCGM", 
          "regime_support_group_maxCGM", "democratic_expericeneCGM", "vote")

stargazer(
  title="Summary Statistics for Electoral Participation Dataset", 
  covariate.labels=c("Male","Age (centered)", "Education (centered)", "Children", "Unemployed",
                     "Income Quintiles","Market Gini (c)","Disp. Gini (c)" , "Clean elections (c)",
                     "Vote Buying (c)", "GDP pc log (c)", "Regime Support Group", "Democratic Experience" ,"Vote"),
  as.data.frame(wvsdata[, cols]), type = "latex", 
  summary.stat = c("n", "mean", "median", "min", "max", "sd")
)

#################################
#Analysis: three-level models: Market Inequality
##################################

v1 <- glmer(vote ~ 1 + (1 | S025) + (1 | S003),
            data = wvsdata, 
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v1)
performance::icc(v1, by_group = TRUE)

### store the intercepts variance, which
### frustratingly, is also called theta in lme4 models
### and it's stored as the sd, hence the need to square it
var_k <- as.numeric(getME(v1, "theta")[2]^2) # level 3 variance
var_j <- as.numeric(getME(v1, "theta")[1]^2) # level 2 variance

### store the alpha value (which lme4 stores as theta = 1/alpha)
alpha <- 1/getME(v1, "theta")

### ICC for level 2
ICC_l2 <- (var_k + var_j)/(var_k + var_j + alpha) 

### ICC for level 3
ICC_l3 <- var_k/(var_k + var_j +alpha ) 



v2 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v2)
anova(v1, v2)
performance::icc(v2, by_group = TRUE)


#macrolevel gini
v3 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + 
              gini_mktCGM + 
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v3)
anova(v2, v3)

#macrolevel gini plus addtional macro controls
v4 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + 
              gini_mktCGM + v2xel_frefairCGM + e_migdppclnCGM + v2elvotbuyCGM +
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
isSingular(v4, tol = 1e-4)
summary(v4)
anova(v2, v4)

# model 4 plus random slopes
v5 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + 
              gini_mktCGM + v2xel_frefairCGM + e_migdppclnCGM + v2elvotbuyCGM +
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 + income_quintilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v5)
anova(v2, v5)

# model 5 plus interactions 
v6 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC +
              gini_mktCGM + v2xel_frefairCGM + e_migdppclnCGM + v2elvotbuyCGM +
              income_quintilesCWC*gini_mktCGM + v2elvotbuyCGM*income_quintilesCWC +
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 + income_quintilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v6)
anova(v5, v6)


# model 6 but with interaction of freedom/fariness and income level and without gini
v7 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC +
              v2xel_frefairCGM + e_migdppclnCGM + v2elvotbuyCGM +
              income_quintilesCWC*v2xel_frefairCGM +  v2elvotbuyCGM*income_quintilesCWC +
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 + income_quintilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v7)
anova(v6, v7)

##Interplot relationship between income and inequality

plot_Interaction_vote_mkt <- interplot(v6, 
                                        var1 = "gini_mktCGM", 
                                        var2 = "income_quintilesCWC", 
                                        ci = 0.95, 
                                        hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of inequality on individual voting") +
  theme_pubr() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "black",
             size = 0.5)
plot_Interaction_vote_mkt
ggsave("output/M1interactioninequality.pdf")


plot_Interaction_vote_buying <- interplot(v6, 
                                        var1 = "v2elvotbuyCGM", 
                                        var2 = "income_quintilesCWC", 
                                        ci = 0.95, 
                                        hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of vote buying on individual voting") +
  theme_pubr() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "black",
             size = 0.5)
plot_Interaction_vote_buying
ggsave("output/M1interactionvotebuying.pdf")


plot_Interaction_vote_cleanelections <- interplot(v7, 
                                          var1 = "v2xel_frefairCGM", 
                                          var2 = "income_quintilesCWC", 
                                          ci = 0.95, 
                                          hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of clean elections on individual voting") +
  theme_pubr() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "black",
             size = 0.5)
plot_Interaction_vote_cleanelections
ggsave("output/M1interactioncleanelections.pdf")

######################################################
#Building Logg Odds ratios of Models


v1_df <- tidy(v1, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v2_df <- tidy(v2, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v3_df <- tidy(v3, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v4_df <- tidy(v4, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v5_df <- tidy(v5, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v6_df <- tidy(v6, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v7_df <- tidy(v7, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")

#Rescaling ordinal and continous predictors by tow standard deviations
v1_df <- v1_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 1") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_quintilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

v2_df <- v2_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 2") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_quintilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))
v3_df <- v3_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 3") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_quintilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

v4_df <- v4_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 4") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_quintilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

v5_df <- v5_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 5") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_quintilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

v6_df <- v6_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 6") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_quintilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

v7_df <- v7_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 7") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_quintilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

#Log Odds Plot Effects Plot

v_models <- rbind(v2_df, v3_df, v4_df, v5_df, v6_df, v7_df)

v_brackets <- list(c("Individual Predictors", "Male", "Income (centered)"), 
                        c("Country-year Predictors", "Market Inequality (centered)", "Regime Support Group (centered)"),
                        c("Cross level Interactions", "Income * Inequality" , "Income * Fairness election"))

plot_main_vote <- {dwplot(v_models,
                          dot_args = list(size = 2, aes(shape = model)),
                          whisker_args = list(size = 0.5, aes(linetype = model)),
                          dodge_size = 1) +
    theme_bw() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("Predicting Voting") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.75, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank()) +
    scale_colour_grey(start = .4, end = .1, # if start and end same value, use same colour for all models 
                      name = "Model", 
                      breaks = c("Model 2", "Model 3", "Model 4", "Model 4", "Model 5", "Model 6", "Model 7"),
                      labels = c("Model 2", "Model 3", "Model 4", "Model 4", "Model 5", "Model 6", "Model 7")) +
    scale_shape_discrete(name = "Model",
                         breaks = c("Model 2", "Model 3", "Model 4", "Model 4", "Model 5", "Model 6", "Model 7"),
                         labels = c("Model 2", "Model 3", "Model 4", "Model 4", "Model 5", "Model 6", "Model 7"))
  
  } %>% 
  add_brackets(v_brackets)

plot_main_vote
ggsave("output/M1dotwhisker.pdf", width = 22 , height = 25 , units = c("cm"))

plot_main_vote <- {dwplot(v_models,
                          dot_args = list(size = 2, aes(shape = model)),
                          whisker_args = list(size = 0.5, aes()),
                          dodge_size = 1) +
    theme_pubr() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("Predicting Voting") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.75, 0.8),
          legend.justification = c(0, 0), 
          legend.title = element_blank())
    
} %>% 
  add_brackets(v_brackets)

plot_main_vote
ggsave("output/M1dotwhisker-color.pdf", width = 22 , height = 25 , units = c("cm"))

#### texreg regression tables main models ####

texreg(list(v1, v2, v3, v4, v5, v6, v7), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)","Male",
                             "Age(centered)","Education(centered)",
                             "Children",
                             "Unemployed",
                             "Income quintiles (centered)",
                             "Inequality (centered)",
                             "Fairness elections (centered)",
                             "GDP pc log (centered)",
                             "Vote Buying (centered)",
                             "Previous Dem. Experience (centered)", 
                             "Regime Support Group (centered)", 
                             "Income * Inequality",
                             "Income * Fairness election",
                             "Income * Vote Buying"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Main Models Voting",
       fontsize = "scriptsize",
       label = "Table:C1", 
       longtable =  TRUE)


###############################
##Using Disposable Inequality##
###############################

## Supplementary Appendix: C1 ##

d1 <- glmer(vote ~ 1 + (1 | S025) + (1 | S003) ,
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d1)

d2 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d2)
anova(d1, d2)

#macrolevel gini
d3 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + 
              gini_dispCGM + 
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d3)
anova(d2, d3)

#macrolevel gini plus addtional macro controls
d4 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + 
              gini_dispCGM + v2xel_frefairCGM + e_migdppclnCGM + v2elvotbuyCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d4)
anova(d2, d4)

# model 4 plus random slopes
d5 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + 
              gini_dispCGM + v2xel_frefairCGM + e_migdppclnCGM + v2elvotbuyCGM +
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 + income_quintilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d5)
anova(d2, d5)

# model 5 plus random slopes
d6 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed+ income_quintilesCWC +
              gini_dispCGM + v2xel_frefairCGM + e_migdppclnCGM + v2elvotbuyCGM + 
              income_quintilesCWC*gini_dispCGM +  income_quintilesCWC*v2elvotbuyCGM +
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 + income_quintilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d6)
anova(d5, d6)

# model 5 without Inequality and interaction
d7 <- glmer(vote ~ sex + ageCWC+ educationCWC + children + 
              unemployed+ income_quintilesCWC +
              v2xel_frefairCGM + e_migdppclnCGM + v2elvotbuyCGM + 
              income_quintilesCWC*v2xel_frefairCGM +  income_quintilesCWC*v2elvotbuyCGM +
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 + income_quintilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d7)
anova(d5, d7)

##Interplot relationship between income and inequality

plot_Interaction_vote_disp <- interplot(d6, 
                                         var1 = "gini_dispCGM", 
                                         var2 = "income_quintilesCWC", 
                                         ci = 0.95, 
                                         hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of disposable inequality on individual voting") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_vote_disp
ggsave("output/M2interactioninequality.pdf")


plot_Interaction_vote_disp_vote_buy <- interplot(d6, 
                                        var1 = "v2elvotbuyCGM", 
                                        var2 = "income_quintilesCWC", 
                                        ci = 0.95, 
                                        hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of vote buying on individual voting") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_vote_disp_vote_buy
ggsave("output/M2interactionvotebuying.pdf")


plot_Interaction_vote_disp_cleanelections <- interplot(d7, 
                                                 var1 = "v2xel_frefairCGM", 
                                                 var2 = "income_quintilesCWC", 
                                                 ci = 0.95, 
                                                 hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of clean elections on individual voting") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_vote_disp_cleanelections
ggsave("output/M2interactioncleanelections.pdf")


#### Plotting Dot-Whisker Plot main Models Vote ####

d1_df <- tidy(d1, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d2_df <- tidy(d2, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d3_df <- tidy(d3, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d4_df <- tidy(d4, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d5_df <- tidy(d5, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d6_df <- tidy(d6, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d7_df <- tidy(d7, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")

#Rescaling ordinal and continous predictors by tow standard deviations
d1_df <- d1_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 1") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_quintilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

d2_df <- d2_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 2") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_quintilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))
d3_df <- d3_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 3") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_quintilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

d4_df <- d4_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 4") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_quintilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

d5_df <- d5_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 5") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_quintilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

d6_df <- d6_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 6") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_quintilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

d7_df <- d7_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 7") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2xel_frefairCGM = "Fairness election (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       v2elvotbuyCGM = "Vote Buying (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2xel_frefairCGM" = "Income * Fairness election",
                       "income_quintilesCWC:v2elvotbuyCGM" = "Income * Vote Buying"))

#Log Odds Plot Effects Plot

v_models <- rbind(d2_df, d3_df, d4_df, d5_df, d6_df, d7_df)

v_brackets <- list(c("Individual Predictors", "Male", "Income (centered)"), 
                   c("Country-year Predictors", "Disp. Inequality (centered)", "Regime Support Group (centered)"),
                   c("Cross level Interactions", "Income * Inequality" , "Income * Fairness election"))

plot_main_vote_disp <- {dwplot(v_models) +
    theme_bw() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("Predicting Voting") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.85, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank())} %>% 
  add_brackets(v_brackets)

plot_main_vote_disp
ggsave("output/M2dotwhisker.pdf", width = 22 , height = 25 , units = c("cm"))


#### Texreg regression table ####

texreg(list(d1, d2, d3, d4, d5, d6, d7), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)","Male",
                             "Age(centered)","Education(centered)",
                             "Children",
                             "Unemployed",
                             "Income quintiles (centered)",
                             "Inequality (centered)",
                             "Fairness elections (centered)",
                             "GDP pc log (centered)",
                             "Vote Buying (centered)",
                             "Previous Dem. Experience (centered)", 
                             "Regime Support Group (centered)", 
                             "Income * Inequality",
                             "Income * Fairness election",
                             "Income * Vote Buying"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Main Models Voting (Disp. Inequality)",
       fontsize = "scriptsize",
       label = "Table:C2", 
       longtable =  TRUE)